function [accuracy, varacc] = filterbank_canonical_correlation_analysis(data,freq,nHarmonics,epochLength,stepTime,Fs,N,a,b)
%FBCCA
%使用方法：放入一种情形下的所有数据，得到正确率
%---input---
%data：数据，size1人，size2刺激频率，size3不同试次，size4通道，size5数据点
%freq：所有刺激频率，单位Hz
%nHarmonics：谐波次数
%stimTime：窗长，单位s
%stepTime：滑窗长度，单位s
%Fs：采样率，单位Hz
%maxFreq：滤波器组频率上界，单位Hz
%N：滤波器组数目
%a：参数指数权重
%b：参数常数项权重
%---output
%accuracy：正确率
%var：关于人的SEM

L = round(epochLength*Fs);
t = (1:L)/Fs;
subjectCount = size(data,1);
stimulusCount = size(data,2);
trialCount = size(data,3);

%% 构造参考信号
for f = 1:length(freq)
    for nh = 1:nHarmonics
        s1 = sin(2*pi*nh*freq(f)*t);
        s2 = cos(2*pi*nh*freq(f)*t);
        refData(:,2*nh-1:2*nh,f) = cat(2,s1',s2');
    end
end

%% 构造滤波器组
[B, A] = CreateFBFilter(Fs);
for n = 1:N    
    w(n) = n^(-a) + b;
end

%% CCA Algorithm
for subject = 1:subjectCount
    for stimulus  = 1:stimulusCount
        for trial = 1:trialCount
            epoch = data(subject,stimulus,trial,:,:);  %单一trial数据
            
            %移除错误数据,根据数据的不同也不同
%             if((subject == 2) && (trial == 10) && (stimulus == 2) && (condition == 2))
%                 epoch = data(session,stimulus,trial-1,:,:);
%             end
            
            epoch = reshape(epoch,size(data,4),size(data,5));
            
            windowCount = floor((size(epoch,2) - epochLength * Fs + 1) / (stepTime * Fs));
            for window = 1:windowCount
                wdata = epoch(:,int32((window-1) * stepTime * Fs + 1):int32((window-1) * stepTime * Fs + epochLength * Fs));
                rho = zeros(N,size(data,2));
                for n = 1:N
                    fdata = filtfilt(B{n},A{n},wdata');
                    for f = 1:length(freq)
                        [~,~,D] = canoncorr(fdata, refData(:,:,f));
                        rho(n,f)=sum(D);
                    end
                end
                rho_sum = sum(w(n)' .* rho.^2,1);
                result(subject,stimulus,trial,window) = find(rho_sum == max(rho_sum));
%                 display(['Subject: ', num2str(subject), 'epochLength: ', num2str(epochLength), ' Stimulus: ', num2str(stimulus),' Trial: ', num2str(trial) ,...
%                     ' Window: ', num2str(window),'Result: ', num2str(result(subject,stimulus,trial,window))]);
            end
        end
    end
end

%% 计算正确率
for subject = 1:size(data,1)
    for stimulus = 1:size(data,2)
        hit = sum(result(subject,stimulus,:,:) == stimulus,'all');
        acc(stimulus) = hit / (size(result,3) * size(result,4)) * 100;
    end
    accuracy(subject) = mean(acc);
end
varacc = sqrt(var(accuracy) / size(data,1));
accuracy = mean(accuracy);
end

